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Abstract 

We propose an efficient Monte Carlo algorithm for simulating a "hardly- 
relaxing" system, in which many replicas with different temperatures are si- 
multaneously simulated and a virtual process exchanging configurations of 
these replica is introduced. This exchange process is expected to let the 
system at low temperatures escape from a local minimum. By using this 
algorithm the three-dimensional itJ Ising spin glass model is studied. The 
ergodicity time in this method is found much smaller than that of the multi- 
canonical method. In particular the time correlation function almost follows 
an exponential decay whose relaxation time is comparable to the ergodicity 
time at low temperatures. It suggests that the system relaxes very rapidly 
through the exchange process even in the low temperature phase. 
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I. INTRODUCTION 



The low temperature phase of spin glasses (SG) and other complex systems generally 
have numerous local minima which are separated to each other by energy barriers. In 
studying such systems, we have to take account of each configuration for these local minima 
and the fluctuation around it. The characteristic time in which the system escapes from a 
local minimum, however, increases rapidly as temperature decreases. This situation causes 
"hardly-relaxing" problem in using conventional Monte Carlo (MC) simulations based on a 
local updating. 

So far, various improvements have been (and is being) made on MC algorithms to over- 
come this problem. There are mainly two categories in such algorithms, i.e., the cluster 
algorithm and the extended ensemble method. The Swendsen-Wang []T| and the Wolff 
algorithms are in the first category in which non-local updating is employed to let the system 
easily change over energy barriers. The second one is based on the idea that an extension 
of the canonical ensemble can be introduced so as to perform MC sampling efficiently at 
low temperatures. The multi-canonical method and the simulated tempering are 
two such algorithms. Some applications of the SG system have been attempted for the 
multi-canonical method [§,^3 and for the simulated tempering [jlO|,|lT|. In this work we 
propose another efficient algorithm belonging to the second category. It can be regarded as 
an parallelized version of the simulated tempering, although the extended ensemble used is 
very different. 

Our method treats a compound system which consists of many replicas of the system 
concerned. The temperature attributed to each replica is distributed in a range including 
both high and low temperature phases. To this system we perform a standard MC simulation 
in the following way: (1) Each replica is simulated simultaneously and independently as 
canonical ensemble for a few MC steps(MCS). (2) Exchange of the configurations for a pair 
of replicas is tried by referring to the energy cost of the whole system. In effect, state 
traverses over the replicas, due to the exchange process, being cooled down and warmed up, 
and thus replica at low temperatures can escape from a local minimum very easily. Since 
the present method never violates the detailed balance in the canonical sense, each replica 
is guaranteed to be equilibrated to the canonical distribution with its own temperature. 
An advantage of our method to the simulated tempering is that we need not to estimate 
weighting factor (or free energy) to equidistribute the probability of visiting temperatures. 

We then apply the exchange MC method to a three-dimensional ± J Ising spin glass(SG) 
model and show that it works well even in the SG phase. To characterize the time scale of 
relaxation in our method, we observe the ergodicity time, which is defined as the average 
MCS per one travel over the whole temperature range, and a modified autocorrelation 
function. It is found that the autocorrelation function follows the exponential decay and 
the relaxation time at low temperatures takes a moderate value which is comparable with 
the ergodicity time. The order parameter distribution P{q) is an even function in the whole 
temperature range, which provide us an evidence for reaching the equilibrium. 

This paper is organized as follows: In Sec. 2 we describe the exchange MC algorithm in 
detail. In Sec. 3 we show how to determine the temperature mesh in order that the exchange 
process occurs properly. In Sec. 4 we report the relaxational behavior of our method by 
applying it to the SG model. The last section is devoted to discussion and summary. 
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II. EXCHANGE MC METHOD 



In our method we treat a compound system which consists of non-interacting M rephcas 
of the system concerned. The m-th rephca, described by a common Hamiltonian 7i(X), is 
associated with inverse temperature jSm, i-e., each rephca is in contact with its own heat 
bath having different temperature (for convenience we take Pm < Pm+i)- A state of this 
extended ensemble is specified by {X} = {Xi,X2, ■ ■ ■ ,Xm}, and the partition function is 
given as 

M M 

Z = Trix} exp(- J] = [] ^iPm), (2.1) 

m=l m=l 

where Z{(3) is the one for the original system. For a set of temperature {/?} given, the 
probability distribution of finding {X} becomes 

M 

P({X,/3}) = n^cq(X^,/3^^), (2.2) 

m 

where 

Peq(X,/3) = Z~\f3)exp{-f3n{X)) (2.3) 

In constructing a Markov process for exchange MC we introduce a transition matrix 
W{X, Prn\X' , Pn) which is a probability of exchanging configurations of the n-th and m-th 
replicas. In order that the system remains at equilibrium, it is sufficient to impose the 
detailed balance condition on the transition matrix: 

P(- ■ ■ ; X, ■ ■ ■ ; X', /5„; ■ ■ ^^^(X, fijX', /5„) 

= P(- ■ ■ ; X', ■ ■ ■ ; X, (3^; ■ ■ ■)WiX', /?^|X, (3^). (2.4) 



exp(-A), (2.5) 



From eq.( |2.3| ) we obtain 

W{X,(3JX',(3,,) 
W{X\l3m\X,f3n) 

where 

A = (/5n - f3n.){n{X) - n{X')). (2.6) 

Therefore the replica-exchange part of transition probability can be expressed as 

if one adopts the Metropolis method. 

For the actual MC procedure, the following two steps are performed alternately: 

1. Each replica is simulated simultaneously and independently as canonical ensemble for 
a few MCS by using a standard MC method. 
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2. Exchange of two configurations and X^n+i, is tried and accepted with the proba- 

bihty W{Xm,/3m\Xm+l:Pm+l)- 

Here we restrict the rephca-exchange to the case n — m + 1 because the acceptance ratio of 
the exchange trial decreases exponentiaUy with the difference (3m — Pn- 

The canonical expectation value of a physical quantity A is measured in usual way as 
follows: 

1 j Vmc s 

= E A{Xm{t)). (2.8) 

-''mcs (=1 

Another expression can be obtained when the exchange procedure mentioned above is re- 
garded as for temperature, i.e.. temperatures, instead of configurations, of a pair of replicas 
are to be exchanged. Then the above quantity is expressed as 

1 Nmcs M 

{A)p ^^Y.T. A{Xmmp,p^(t), (2.9) 

-'Vmcs t=\ m=l 

where we introduce the time-dependent inverse temperature (3m{t) and the configuration Xm 
in this temperature-exchange scheme. Note that both schemes are completely equivalent to 
one another. One can choose either of the two schemes in actual implementation of the 
present method. 



III. DETERMINATION OF TEMPERATURE 

The parameters we have to determine are only the set of (inverse) temperatures {/5m}- 
The highest temperature should be set in the high temperature phase where relaxation 
(correlation) time is expected to be very short and there exists only one minimum in the 
free energy space, otherwise the system would not completely forget where it was trapped 
before even if it visits to the highest temperature. On the other hand, the lowest temperature 
is somewhere in the low temperature phase whose properties we are interested in. In this 
sense the temperature range is considered to be given. Then we estimate the number of 
temperatures required in the range as follows. In order that each replica wanders over the 
whole temperature region the acceptance probabilities of the exchange process for every pair 
of rephcas at different temperatures have to be of order of one and nearly constant. The 
logarithm of the probabihty e~^ of exchanging /3„ and Pn+i — Pn + ^ ^sto order 6'^ 

A = S{H{Xn+i) - H{Xn)) - S^-^E, (3.1) 

where the instantaneous value of the energy Ti is approximated by the thermal expectation 
value E. Since the energy E is an extensive variable, S should be of order of to satisfy 
the condition that A ~ 0{1). In other words, the number of temperatures in the range we 
have to simulate is of the order of V^. In the case where a second-order phase transition 
exists, the number is modified as viV^^, where a and p are the exponent of the specific 
heat and the correlation length. This means that we need more replicas if the specific heat 
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diverges at the transition temperature. In spin glass and other glassy systems, the specific 
heat has usually no singularity at the transition temperature. 

A set of temperatures {Pm} can be obtained by the iteration procedure which is first 
introduced to the simulated tempering. |]10[ For given {Pm}, the acceptance ratios {pm} are 



evaluated by simulating an appropriate MCS. Then a new set 1/3^} is constructed by using 
the old set {Pm}'- 

P'l = Pi, 

P'm = P'm-1 + iPm ' Pm-l) — , (m = 1, 2, ■ ■ ■ , M) , 

C 



M 



M 



— E Pm- (3.2) 



m=l 



In each iteration the MCS required to estimate pm are found not so long. We will show a 
concrete examination in the next section. 

In closing this section, we propose some necessary conditions to check the efficiency and 
the equilibration. 

(i) The exchange happens with a non-negligible probability for all adjacent pairs of 
replicas. 

(ii) Each replica moves around the whole temperature range in suitable MCS. 

(iii) In moving the temperature space the system forgets where it was trapped. 

We can obtain the temperatures {Pm} satisfying the condition (i) by the iteration procedure 
mentioned above. The ergodicity time defined in the following section gives a criterion 
for the condition (ii). Unfortunately, even if the condition (i) and (ii) are satisfied, it is 
possible that the condition (iii) is violated. One trivial possibility for this is that the largest 
temperature is not high enough. If some of these conditions are broken down, we may not be 
able to improve the present method due to the absence of any other controlling parameters. 



IV. THE MODEL AND SIMULATIONS 

We consider the three dimensional ±J Ising SG model on the simple cubic lattice. The 
real-replica Hamiltonian is defined as 

n{{a,r}) = - ^ J,^.(a,or, +r,r,), (4.1) 

where {a} and {r} take the values ±1. The interactions {Jij} are quenched random variables 
taking ±1 with equal probability. To be more accurate, {Jij} are distributed so that the 
number of the ferromagnetic bonds, Jij = +1, is exactly a half of the total bonds. Then 



the Edwards- Anderson SG order parameter |jT2[ is computed as the overlap between the two 
copies, 

5 = ^[E(^^^0], (4.2) 



where (. . .) and [. . .] denote thermal average for the Hamiltonian (|4.1|) and random average, 
respectively. We have simulated on lattices of linear size L = 6, 8, 12, and 16 with periodic 
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boundary condition. For local updating we have adopted the conventional sublattice-flipping 
with heat bath method. For the exchange process the replica pairs {Pm, Pm+i) are divided 
into two subgroups, i.e., odd-m and even-m groups. In actual updating the exchange trial 
is performed for one of these subgroups after each local updating. It means that our one 
MC step consists of one local updating MC step and a half exchange-trial per rephca-pair. 



A. temperature setting and the ergodicity time 

The actual local updating have been performed with the multi-spin coding technique 
15Hl7| which simulates 32 different physical systems at once. This fixes the total number 



of temperatures to be M = 32 in all lattice sizes. To determine the temperatures {Pm} a 
typical sample of size L = 12 has been used. The initial condition we used is 

TTl — 1 

Pm=f3l + {f3M-f3l)j^, (4.3) 

with f3i = 0.4 and (3m = 2. We then performed 400 MC steps of simulation for each iteration 
of eq.( p.2|) to evaluate the acceptance probabilities {pm}, and found that the temperatures 
converge after several iterations. Finally we smoothed by spline interpolation and 

obtained the temperature range (0.86 < T/J < 2.39) which includes both the high and low 
temperature phases. These temperatures so obtained was used for all samples and sizes. In 
fig. we show the acceptance ratios for various sizes. These values are found to be of order 
of one and do not depend on bond realization so sensitively. 

In order to confirm that the obtained temperature set satisfies the condition (ii) men- 
tioned in the previous section, we investigated the ergodicity time te which is defined as the 
average MC step for a specific replica to move from the lowest to the highest temperature. 
The dependence of our observed ergodicity time on the lattice size is shown in Fig. |^, with 
that of the multi-canonical method by Berg et al. If the total number of temperature 
points we simulate is fixed as in the present analysis, the acceptance ratio, e~^, of exchange 
process between replicas reduces exponentially as a function of system size L and the te 
grows rapidly with increasing L. As shown in Fig. ^, however, it is clear that our ergodicity 
time is much shorter than the multi-canonical ergodicity time in a lattice size suitable for 
our actual simulation. 



B. relaxation 

Because each replica wanders over the temperature space, the (equilibrium) time cor- 
relation function can not be argued in the ordinary sense. Here to investigate relaxation 
dynamics of the present method, we study an autocorrelation function in the temperature- 
exchange scheme, 

?(^,/5J = ^E(^^(0)^^W)pl, (4.4) 

i 

where (■ ■ ■)p™th nieans an average along trajectories in temperature space evolving from 
the initial state {crj(O)} with temperature Pm- The function q(t,P) is expected to include 
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the slowest relaxation mode of the present method. In Fig. ^ we show q{t, j3) at various 
temperatures for L = 12. The data are averaged over 10 samples. As indicated by the 
straight lines, the function q{t,(3) follows nearly exponential decay for large time t. We 
evaluate the largest relaxation time from least square fit to a single exponential 

~goexp(-t/r) t > 1. (4.5) 

As shown in Fig. ^(a) the relaxation time r has a crossover temperature ~ 1.5, above 
which r behaves like the relaxation time obtained by a conventional local spin flip dynamics. 
Below Tcr, on the other hand, r seems to almost become saturated and to be comparable with 
te, and the relaxation amplitude go grows rapidly (see Fig. §(b)). It suggests that the largest 
relaxation mode of the exchange dynamics dominates the relaxation at low temperatures. In 
any case we can see that the system relaxes within a reasonable MCS and thus the condition 
(iii) is satisfied. 

C. order parameter distribution 

In this subsection we show the obtained distribution function Pj{q) of overlap function 
for a sample, which is defined by 

PAQ) = {Siq-j^j:^-^n))- (4-6) 

i 

The distribution function Pj{q) is important to SG study because all physical quantities of 
our interest can be obtained from this function. The distribution is expected to be an even 
function from the invariance of the Hamiltonian ( |4.1| ) under the transformation cXi — o"j. 
We demonstrate in Fig. |^ that for L = 16 the distribution Pj{q) is symmetric even below 
the SG transition temperature. Note that only 3 x 10^ MCS, which is about 10 times larger 
than Te are used to obtain Fig. |. The distribution of the same sample obtained by the 
conventional MC method is also shown in Fig. |^ with the broken line. Obviously, it is far 
from symmetric and even the peak positions differ with each other. It strongly suggests that 
the system does not reach the thermal equilibrium yet by the conventional MC method. 

V. DISCUSSION AND SUMMARY 

One may consider that, for random systems, should be determined for each sample 

separately. It is, however, not the case because the exchange probability {pm} is insensitive 
to bond realization as long as Pm+i — is sufficiently small. It is contrasted to the case of 
the multi-canonical method, where the energy density should be estimated for each sample 
because it strongly depends on the structure of local minima in a sample. |^ 

An advantage of the present method to the simulated tempering is that we need not 
estimate weighting factor (or free energy) to equidistribute the probability of visiting tem- 
peratures. This is because the phase space of the present method is a direct product of the 
original one in contrast with that of the simulated tempering, which is given by a direct 
sum. Since the weight factor in the simulated tempering is an extensive parameter, the 
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probability is very sensitive to it, and temperature traversing would be easily broken if the 
estimation is not so good. 

A more important merit of the present method appears in evaluating P{q) of the SG sys- 
tem, where two copies of a system are independently simulated. We have M replicas for each 
copy, so that we can take M{M — l)/2 un-correlated samplings of overlap in one simulation 
as compared with M sampling is by M independent runs of the simulated tempering. 

As shown in the previous section, the model to be simulated is not supposed to have any 
specific aspect, so that we can, in principle, apply the present method to various models 
without modification. Unfortunately the present method seems not to work well in the 
systems exhibiting the first order phase transition. Since the energy has a finite gap at 
transition point, the exchange between replicas below and above transition temperature 
hardly occurs in a large lattice, meaning that the condition (i) in section 3 is not satisfied. 

In summary we have proposed a new MC method for simulating hardly relaxing systems. 
We have applied this method to the three dimensional ± J Ising spin glass system, and found 
that the system really traverses over wide temperature space and the largest relaxation 
time in this dynamics is given by the ergodicity time, which is much smaller than the 
conventional one. As a result, the order parameter distribution P{q) can be obtained below 
Tc up to L = 16 within much shorter than the conventional MCS. Numerical results of 
physical quantities and discussion on the SG phase transions will be discussed in detail 



elsewhere. \2G 
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FIGURES 



FIG. 1. Temperature dependence of the acceptance probability of exchange process. The sam- 
ple-to-sample errors are smaller than 1/10 of the size of symbols. 

FIG. 2. Ergodicity time te as a function of the system size V = L^. The circles correspond to 
our data, and squares to multi-canonical method by Berg et al. The unit of relaxation time is MC 
step per spins. 

FIG. 3. The modified autocorrelation function q{t,l3) for 1/(3 = T = 0.862, 1.088, 1.256, 1.399, 
and 1.609 ( top to down). 

FIG. 4. Temperature dependence of the relaxation time r (a) and its amplitude qo (b) of the 
modified autocorrelation function (Fig. 

FIG. 5. Distribution function Pj{q) of the overlap for a sample with L = 16 at T = 0.924. 
The solid and dashed line represent P{q) by using the exchange MC and by a conventional MC 
method, respectively. The same MCS are used in two methods. 
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